function phi = phi_3D(x,y,z,h,M,r,theta)
%%compute the potential of a cylindrical PM
%%[x,y,z]:interested point
%%[r,theta]:integral variable

phi_hat = (r.*M)./sqrt((r.*cos(theta)-x).^2+(r.*sin(theta)-y).^2+(h/2-z).^2)...
         -(r.*M)./sqrt((r.*cos(theta)-x).^2+(r.*sin(theta)-y).^2+(-h/2-z).^2);

phi = phi_hat/(4*pi);
